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Interactive multiple anisotropic scattering in clouds 
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^EVASION - LJK / Grenoble Universites - INRIA Davis 



Abstract 

We propose an algorithm for the real time realistic simulation of 
multiple anisotropic scattering of light in a volume. Contrary to 
previous real-time methods we account for all kinds of light paths 
through the medium and preserve their anisotropic behavior. 

Our approach consists of estimating the energy transport from the 
illuminated cloud surface to the rendered cloud pixel for each sep¬ 
arate order of multiple scattering. We represent the distribution of 
light paths reaching a given viewed cloud pixel with the mean and 
standard deviation of their entry points on the lit surface, which we 
call the collector area. At rendering time for each pixel we deter¬ 
mine the collector area on the lit cloud surface for different sets of 
scattering orders, then we infer the associated light transport. The 
fast computation of the collector area and light transport is possi¬ 
ble thanks to a preliminary analysis of multiple scattering in plane- 
parallel slabs and does not require slicing or marching through the 
volume. 

Rendering is done efficiently in a shader on the GPU, relying on 
a cloud surface mesh augmented with a Hypertexture to enrich the 
shape and silhouette. We demonstrate our model with the interac¬ 
tive rendering of detailed animated cumulus and cloudy sky at 2-10 
frames per second. 

1 Introduction 

The realistic rendering of clouds is still a challenging problem. 
Multiple anisotropic scattering must be simulated in a volume, 
the lack of absorption makes the convergence very slow, and de¬ 
tailed clouds require high resolution volumes. In the scope of 
real-time rendering this is even worse: volume radiosity or Monte 
Carlo methods cannot converge in real-time, and volume render¬ 
ing cannot be achieved with enough resolution for sharp clouds 
such as cumulus. Precomputations forbid the animation of clouds 
and light source, and most real-time models do not account for 



view-dependent effects or miss important visual features such as 
backscattering. 

To solve these problems, our approach represents clouds as 
surface-bounded volumes and optimizes the calculation of light 
transport inside the cloud from the illuminated cloud surface to the 
rendered cloud pixels. For this, we study and characterize the light 
transport for each order of scattering. 

Our contributions are: 

- a new model characterizing the light transport (i.e., amount of 
transmitted energy) between the fiat surface of a slab and a 
point p anywhere in the slab at any given order of scattering 
(see Figure 1(a)); 

- a new model characterizing the distribution on a slab surface of 
the entry points of the light rays reaching any given point p in 
the slab for any given order of scattering. We call this entry area 
the collector area (see Figure 1(a)) ; 

- an iterative algorithm to determine this collector area on an arbi¬ 
trary cloud surface for any given rendered point p, and to com¬ 
pute the associated light transport (see Figure 1(b)) ; 

- An efficient representation of detailed cloud shape using a Hy¬ 
pertexture [Perlin and Hoffert 1989] on a surface mesh (see Fig¬ 
ure 1(c)) ; 

- A GPU implementation of these contributions, resulting in 
highly detailed rendering of animatable clouds interactively 
with all-orders multiple Mie scattering. 



(a) (b) (c) 

Figure 1: Overview of our contributions, (a): In an preliminary analysis, 
we characterize light transport in a slab via a collector area representing 
location of incoming light for each scattering order, (b): We use this char¬ 
acterization to find collectors on an arbitrary cloud shape and compute the 
light transport, (c): Our cloud model is represented by a Hypertexture, with 
procedural details on the boundary and a homogeneous core. 

Contrary to most previous real-time models our multiple scattering 
simulation: 

- does consider the high orders of scattering which are responsi¬ 
ble for diffusion and backscatter, as well as the low orders of 











scattering which are responsible for the glory, the silver lining 
on the silhouette, and the appearance of thin parts; 

- uses the physieally based strongly anisotropie Mie phase func¬ 
tion (not Rayleigh, Gaussian or isotropie) ; 

- does not need to walk through the volume ; 

- does not rely on any shape-dependent preeomputation. 

Representing the cloud shape with a Hypertexture on a surfaee 
mesh allows us to effieiently render inhomogeneous boundaries as 
well as sharp detailed clouds such as cumulus, contrary to meth¬ 
ods based on sliced volumes. We can thus render sharp, fluffy or 
complex wispy clouds at little cost. 



Figure 2: Log plots (inset: polar log plots) of commonly used phase func¬ 
tions. Red: Rayleigh. Green: Henyey-Greenstein with g = .99. Blue: Mie. 

2 Cloud physics 


2.1 Density and size of droplets in clouds 

Real convective elouds are not blurry on boundaries, they are usu¬ 
ally sharp or wispy (see Figures 5(e), 5(d)). Collapsing eloud tur¬ 
rets or weak clouds may have a larger wispy layer. These inhomo¬ 
geneities in the eloud liquid water content {i.e., mass density) are 
a strong visual feature of clouds. Inside the cloud, this density can 
also vary - especially near rain condition - due to eoalescenee of 
droplets. 

The size of droplets is charaeterized at any loeation by a droplet size 
distribution (DSD), which is generally modeled in literature by a 
lognormal or a modified Gamma funetion [Levin 1958]. The optical 
properties of a cloud depend highly on the droplet size. Accounting 
for the DSD radieally ehanges the resulting phase funetion, thus the 
visual appearance. 

Still, these data within a cloud is not often available and its physies 
{e.g., the eoalescenee mechanism or the evolution of the DSD) is not 
fully understood. Since computer graphics applications require pri¬ 
marily visual plausibility, approximations of these values are com¬ 
mon. As an example, density variations are visually more impor¬ 
tant on the clouds boundaries -where they are directly visible- than 
in the eloud core -where they only influence the appearance indi¬ 
rectly, through high-order multiple scattering. Usually, the DSD 
over a cloud (if using a DSD at all) is assumed eonstant. 

2.2 Phase function 

The phase function of a cloud droplet is given by the Mie theory. 
Common approximations are Gaussian or Henyey-Greenstein func¬ 
tions. As seen in Figure 2 the Mie funetion combines a strong nar¬ 
row forward peak (51% of energy), a wide forward lobe (48% of 
energy), and a eomplex baekward lobe with peaks. These three fea¬ 
tures all yield very specific and visible effects at the seale of the 



Figure 3: Some results of our light transport analysis. Left: BSDF for 
a point of view at a depth d = 100m inside a cloud slab of thickness t = 
1000m, with an incident illumination angle 4 >l = 75°. Areas represent the 
contribution of different orders of scattering. Even high orders (up to 30 
scattering events) show an anisotropic behavior, while orders > 30 show 
an isotropic behavior. Right: Contribution of different orders of scattering 
for the reflectance (in red) and transmittance (in blue) of a slab of varying 
thickness. Isotropic orders (31-oo) begin to appear at thicknesses > 100m. 
Anisotropic orders (1-30) play a role up to 1000m in transmittance and 
contribute to 30% — 100% of the reflectance. 

cloud. The absenee of baekward peaks means no glory or fogbow. 
No narrow forward peak means huge underestimation of global 
transmittance and anisotropy. Gaussian and Henyey-Greenstein 
functions do encode a lobe and ease ealculation but they give far 
from aceurate visual effects. These approximations miss eloud fea¬ 
tures sueh as the glory and fogbow, and blur out the narrow forward 
scattering peak. Rayleigh scattering is even less appropriate since 
it is symmetrical (50% backward) and it physically corresponds to 
molecular scattering (giving its blue color to the sky), as opposed 
to scattering due to droplets. 

In the visible speetrum, water cloud albedo can be eonsidered as 1 
sinee there is no absorption (all light is either refieeted or transmit¬ 
ted). Note that some atmospheric phenomena commonly credited 
to elouds are aetually eaused by other elements (e.g., rainbows are 
eaused by rain, and sundogs by ice crystals in the atmosphere). 

2.3 Anisotropy 

Multiple scattering is strong in clouds and often shows anisotropic 
behavior. Clouds span hundreds to thousands of meters and the 
mean free path of a light ray is about 20m. As a result, most rays 
will be seattered multiple times before exiting the cloud. Because 
of the highly anisotropic nature of the Mie phase function (99% of 
the light is seattered in the forward direction), even multiple scat¬ 
tering can be anisotropic. From our analysis of light transport in a 
slab, we estimate that the light behavior is isotropie only after about 
30 scattering events. By isotropic, we mean that multiple scatter¬ 
ing behaves as if the phase funetion of the medium was isotropic, 
not that the resulting Bidirectional Scattering Distribution Function 
(BSDF) is isotropic. 

3 Previous Work 

Simulating multiple anisotropic Mie seattering through Monte- 
Carlo integration is not practieal in real-time. Previous optimiza¬ 
tion approaches rely on various simplifieations, whieh we discuss 
below. 

3.1 Phase function 

Sinee the Mie phase function is complex and expensive to com¬ 
pute, it is often approximated using other funetions sueh as Henyey- 
Greenstein [Max 1994], Gaussian [Premoze et al. 2004] or even 
Rayleigh [Harris and Lastra 2001] (see comments in Section 2.2). 
[Riley et al. 2004; Bouthors et al. 2006] preeompute the Mie 



























Figure 4: Various kinds of light paths in participating media, brought by 
different orders of scattering. Low orders of scattering bring short, steep¬ 
turning paths (red). Higher orders bring long, slow-turning paths (green). 
Highest orders bring complex, spread, diffusive paths (blue). The Most 
Probable Paths approach [Premoze et al. 2003] reproduces the green ones, 
but underestimates red and blue ones. 

phase function for a given DSD, which reproduces real-life fea¬ 
tures (glory, fogbow, etc.). We use the same approach for our 
phase function. 

3.2 Light transport 

[Kajiya and von Herzen 1984] and [Blinn 1982] consider either 
low albedo or low density: only single scattering is considered. 
This neglects all multiple scattering effects. Assuming the trans¬ 
port is mostly forward allows for real-time single-pass algorithms 
such as slice accumulation [Dobashi et al. 2000; Harris and Lastra 
2001; Riley et al. 2004] but this neglects some of the multiple scat¬ 
tering effects such as backscattering. Computing multiple scatter¬ 
ing [Nishita et al. 1996] for only the lower orders only has the same 
effect. The diffusion approximation for multiple scattering [Stam 
1995; Jensen et al. 2001] allows efficient computations of high or¬ 
ders but neglects anisotropy in multiple scattering. 

[Premoze et al. 2003] introduced the idea of Most Probable Paths 
(MPP) in participating media. The root idea is that most of the 
photons arriving at one point in one direction roughly followed the 
same path. Thus, integrating light transport only along this path 
is sufficient to account for most of the energy. In addition, [Pre¬ 
moze et al. 2004] speed up this technique and account for the spatial 
spreading of light around this mean path through an analytical for¬ 
mulation. This approach has been brought to real-time [Hegeman 
et al. 2005] using graphics hardware for volume slicing in a man¬ 
ner similar to [Harris and Lastra 2001; Riley et al. 2004]. These 
methods based on two slicing passes (accumulating the fiux from 
the light source in the slices, then from the slices to the eye) restrict 
the variety of light paths accounted for. 

As mentioned in [Premoze et al. 2003; Hegeman et al. 2005] the 
main limitation of the MPP approach is that the paths it computes 
are mainly of high order, with a small angle change per scatter¬ 
ing event. Thus, paths of low order as well as diffusive paths are 
underestimated. These paths contribute to most of the lighting in 
thin cloud parts and in the backscattering, and thus should not be 
neglected (see Figure 4). 

We address these limitations by treating light paths of all orders. 
We also propose a new, faster light transport computation approach 
that does not rely on slicing or marching through a volume. Our 
approach is also inspired by the radiative transfer studies on ba¬ 
sic shapes such as slabs [Chandrasekhar 1960; Krim and El Wakil 
1986; Mobley 1989; Max et al. 1997]. We use hardware-friendly 
representations such as depth maps [Dachsbacher and Stamminger 
2003] to implement our rendering algorithm on the GPU. 



(a) [Harris and Lastra 2001] (b) [Riley et al. 2004] 



(c) Photograph (d) Photograph 



(e) [Schpok et al. 2003] (f) [Gardner 1985] 


Figure 5: Top: real-time (a) and interactive (b) CG cumulus clouds using 
billboard or slices. Middle: wispy (c) and sharp (d) real cumulus clouds. 
Bottom: CG clouds using 3D noise (e) and surfaces with procedural de¬ 
tails (f). Adding procedural details gives a less blurry and more contrasted 
appearance, increasing realism. 

3.3 Cloud densities representation 

Since real cloud density data is difficult to measure or simulate, 
early work used procedural models [Gardner 1985; Dobashi et al. 
2000]. Recent approaches have been using fluid simulation tech¬ 
niques [Harris and Lastra 2001] or atmospheric data [Trembilski 
and BroBler 2002; Riley et al. 2004]. However, fluid simulation 
can only be computed at coarse level for real-time applications, 
and atmospheric data have low resolution. As a result, adding 
high frequency details is necessary to avoid a blurry appearance 
(see Figure 5). [Ebert 1997] combines Perlin solid noise [Perlin 
1985] and implicit surfaces. [Schpok et al. 2003] advect a noise 
texture [Neyret 2003], but do not account for the noise in the light 
transport computation, which gives them a uniform appearance (see 
Figure 5(e)). These ideas inspired us to combine two representa¬ 
tions -meshes and 3D textures- to model the clouds at two different 
scales. 

Various rendering primitives have been considered to render clouds. 
Since clouds are a 3D distribution of droplets, volume grids have 
often been used to describe them and volume rendering techniques 
to render them. Methods using billboards or volume slices [Harris 
and Lastra 2001; Premoze et al. 2004; Riley et al. 2004; Dobashi 
et al. 2000] result in a lot of overdraw^ which is computationally 
expensive. While textured slices are an efficient way of rendering 
gaseous phenomena, they are not the best option for dense clouds. 


^i.e., a given pixel is rasterized many time by different rendered primi¬ 
tives. 














(b) 



Figure 6: Our cloud representation. A mesh is used to describe the outer 
cloud boundaries at low resolution. A procedural volumetric hypertexture 
adds details under the boundary up to a certain depth h inside the cloud. 
The core is considered homogeneous. 

where most pixels of the back slices are hidden by those in front and 
therefore wasteful to render. Moreover, 3D texture memory limits 
the resolution of such models, resulting in blurry silhouettes and a 
lack of details (see Figures 5(a), 5(b)). 

Since cumulus clouds are dense and often have a sharp interface, 
they also have been represented as surface-bounded volumes such 
as sets of ellipsoids [Gardner 1985; Elinas and Stiirzlinger 2000] or 
meshes [Trembilski and BroBler 2002]. To avoid the “hard” appear¬ 
ance of polygonal surfaces, they use a procedural shader simulat¬ 
ing the detailed silhouette and giving a volumetric impression (see 
Figure 5(f)). Light transport through the volume is not simulated. 
[Bouthors et al. 2006] do simulate light transport inside a mesh, but 
they consider no enrichment at all on silhouettes which thus appear 
polygonal and opaque. 

We take advantage of both volumes and surfaces by representing 
the high-scale cloud boundary with a mesh and the high-frequency 
density variations at cloud borders with a Hypertexture [Berlin and 
Hoffert 1989]. The rest of the cloud interior is considered homoge¬ 
neous, so only a thin layer of voxels is necessary (see Figure 6). 

4 Overview of our method 

Our shading approach is based on the analysis of light transport 
for each scattering order. Observing that paths of different orders 
have different anisotropy and spreading behaviors, we treat them 
separately: 

- order 1 (single scattering), which is the most anisotropic and 
depends on the finest details, is computed using an analytical 
form of the scattering equation (see Section 6.2) ; 

- orders 2-oo (multiple scattering) are computed in 8 separate sets 
(2, 3-4, 5-6, 7-8, 9-12, 13-18, 19-30, 31-oo) using our collector- 
based algorithm (see Section 6.1) ; 

- opacity is computed by integrating the extinction function 
through the cloud volume (see Section 6.3). 

Our collector-based approach follows and extends the idea of Most 
Probable Paths (MPP) [Premoze et al. 2003]. We consider one 
most probable path and spreading per set of scattering orders. More 
specifically, we consider a collector area, which is the piece of sur¬ 
face through which enters 95% of the light that reaches the current 
rendered pixel in the view direction (see Figure 7(b)). This col¬ 
lector is defined by its center c on the cloud surface and width a. 
For a given pixel, for each set of scattering orders, we look for this 
collector, and compute the corresponding light transport. 

To find this collector area on the lit cloud surface, we rely on an 
algorithm that iteratively matches the cloud surface with the most 
probable collector area. This algorithm is described in Section 6.1. 






(a) 


Figure 7: (a): Setup and notation for the canonical light transport T 
through a collector area (c, cr) between a lit slab surface and a point p in 
the slab, depending on the input parameters (fvi for each 

set of scattering orders. Note that'ipL is always zero (the reference frame is 
aligned with the light direction), (b): Light transport in the cloud between 
its lit surface and a point p. For each set of scattering orders the incom¬ 
ing light at p is assumed to come from a collector area. We characterize 
this transport by looking for a fitting slab, then relying on the canonical 
transport function. 

By definition, the role of the cloud surface outside the collector is 
negligible, thus we can locally approximate the cloud shape as a 
plane-parallel slab aligned on the collector in order to simplify the 
computation of light transport (see Figure 7(b)). 

Computing anisotropic light transport even for a simple shape like a 
slab is still very complex [Chandrasekhar I960]. To accelerate the 
light transport computation, we characterize the radiative transfer 
in a slab through a canonical transport function. This is described 
in Section 5. We obtain it by analyzing numerous Monte-Carlo 
simulations for different slab parameters (see Appendix A). 

We implement our whole rendering approach on the GPU, using 
depth maps to represent the cloud surface. This is described in Sec¬ 
tion 7. In Section 8 we introduce our enrichment of the cloud model 
and especially of its silhouette in a Hypertexture upon the surface. 
The gathering of all steps within a single shader is synthesized in 
Section 8.1. Then we present our results and performance in Sec¬ 
tion 9. 

5 Characterizing canonicai iight transport 

5.1 Simulation 

In this section we describe our experimental setup in the canon¬ 
ical case of a plane-parallel slab of thickness t (see Figure 7(a)). 
For each set of scattering orders we compute the light transport T 
(i.e., the proportion of transmitted energy), the collector center c 
and the collector standard deviation a. 

These values (T, c, a) are computed against 5 input parameters: 
slab thickness t, viewpoint depth d, viewing angles 
and lighting elevation angle fiL, in the reference frame shown 
in Figure 7(a). We define pv = cos(0v^), pL = cos(0l), 
LJv = view direction, ujl — light direction, cos^ = uy - do l, 
p = (0, —(7, 0) = viewpoint location. Since we want to character¬ 
ize light transport up to any point p within the volume, we consider 
values of the viewpoint depth d within the slab (0 < d < t). View¬ 
points outside the slab correspond io d = 0 or d = t. 

We ran numerous Monte-Carlo simulation of light paths for vari¬ 
ous values of these 5 parameters (see Appendix A for the details). 
We store the results (T, c, a){(t)v , (j)L^d,t) in a 5D table for 

each set of scattering orders. We call (T, c, a){(l)v, fjy, fiL,d, t) 
the canonical transport function. 
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Figure 8: Some resulting BSDF of our light transport analysis. T is plotted against fy ('ipy fixed at 0°) for various input parameters (cpL = 25° j. 
Blue: Monte-Carlo simulations. Red: our fitting. 


5.2 Data compression 


6 Estimating iight transport in ciouds 


These 5D tables (1 per set of scattering orders) encode the canon¬ 
ical transport function, that is, the macroscopic light behavior in a 
slab of cloud depending on the input parameters {<fiy , fjy 
Such a set of 5D tables cannot be stored in GPU memory. We 
compress them by approximating the computed results with empir¬ 
ical functions. Appendix A describes our experimental and fitting 
setup. 

• For the light transport T, this fitting yields 




P-A-X-pl 


log{d + D) - 
\og{B + Df 


(d-sA 

2C‘^ 


with A = A(t, py), B = Bi (t, py) — ^2 {t, Pl), C = c(t, py), 
D = T>(t,py), X — P — P(^), Le., a 5D table for 

T reduces to an analytical function and seven 2D tables. P en¬ 
codes the anisotropy of the result and equals 1 if the light behav¬ 
ior is isotropic (orders 31-oo). X modulates the result accord¬ 
ing to the lighting angle and is inspired from Chandrasekhar’s X- 
function [Chandrasekhar I960]. A, Bi, B 2 , C, D are the param¬ 
eters of a “skewed” Gaussian function encoding the light behavior 
according to the depth of the viewpoint. Figure 8 shows the result 
of this fitting. 

• For the collector center c, our compression results in 

c = (cx,0,Cz) with 

Cx = Ax\og{lE ■ d)Bx 
Cz = Az\og{lE ■ d)Bz 
Ax = Fsin'^v sin(G • 

Bx = Hsmfi)ysin((l)L) 

Az = / + J[cos'0\/sin(Ff • 

Bz = M N cos 'fiy 


with E = E = ¥{py), G = G{py), H = Ji{py), 

I = l{py), J = J(py), K = K(py), L = L{py), 
M = 'M.(py), N = lSi(py, Pl)- 

• Our compression of the collector size a gives 


cr = O + Q • t log(l + R • (i) + S log(l + T • t) 


where O, Q, R, S, T are constant for a given set of scattering 
orders. 


6.1 Multiple scattering 

Considering a given cloud pixel to be rendered (corresponding to a 
location p in the cloud), for each scattering order we look for the 
collector area (c, a), i.e., the origin of dispersed light paths on the 
lit cloud surface which ended at p in eye direction (see Figure 7(b)). 
Since a cloud is a volume of inhomogeneous densities, defining its 
surface is difficult. We choose to define it as the boundary where the 
density p crosses a user-defined threshold value po- As explained in 
the overview, we assume that the cloud behaves locally like a slab 
tangent to the surface at the collector location. Here, surface ori¬ 
entation is a scale-dependent notion. Since we are interested in the 
surface of incoming dispersed light paths, our local surface orienta¬ 
tion is obtained through filtering the cloud surface according to the 
dispersion standard deviation a (dashed cloud shape in Figures 7(b) 
and 9). Sections 7 and 8 give more details on this filtering. 

To find this collector location, we iterate as shown on Figure 9. We 
start at step 1 with a first collector of size cro at location co which we 
project in step 2 along the light direction on the cloud surface at Cq. 
The corresponding slab is tangent to the surface (filtered according 
to (Jo) at Co . View and light parameters (j^L^dA) accord¬ 

ing to this slab are obtained by simple geometric transformations. 
In step 3, the canonical transport function provides us with the col¬ 
lector location ci(0a/, fjy, fih^dA) and size cn^fiy, fjy, fh^dA) 
corresponding to such a configuration, which is likely to be differ¬ 
ent from Cq. We project ci on the cloud surface filtered by cri at 
c'l in step 4, and we iterate up to convergence, i.e., Cn ~ c^. We 
then take (c, d) = (cn, cr^). Once we found the collector, we can 
obtain the corresponding light transport T {i.e., the amount of en¬ 
ergy transmitted from this collector to the eye) with the same input 
parameters. We run this algorithm for each set of scattering orders 
on the GPU (see Section 7). 

This iterative algorithm is in fact a fixed-point method applied on 
c, i.e., we are looking for value x that satisfies f{x) = x, where x 
here is our collector parameters (c, a) and / represents our canoni¬ 
cal transport function. Since the search space of this collector is re¬ 
stricted to the cloud lit surface, this algorithm cannot diverge. How¬ 
ever, like any basic fixed point method, it can reach a state where it 
loops between two values without converging {i.e., Ci A ^i-i 
d = Ci- 2 ). It can also take too large steps and miss the solution. 
To avoid these cases, we limit the size of each step (we constraint Ci 
such that ||ci — Ci_i II < a A. We start the algorithm with an initial 
position Co = p. We use a large initial size cro. Indeed, if the initial 
collector spans the whole cloud, the projected result at step 2 will 
be a first-order approximation of the whole lit cloud surface, which 
























Figure 9: Summary of our iterative algorithm to locate the collector (blue 
area). Solid cloud shape: unfiltered cloud surface. Dashed cloud shape: 
surface filtered according to cr. Dashed box: slab matching the filtered sur¬ 
face at collector position c. Orange boxes represent the use of our canonical 
transport function. 

is a good starting point for our algorithm in terms of convergence 
speed. Subsequent steps will then “refine” the lit cloud surface at 
the most probable collector location. 

6.2 Single Scattering 

For the first scattering order, which constitutes a degenerate case 
for our collector-based algorithm, we simply integrate analytically 
the Mie single scattering along the eye direction. This allows us to 
account for the inhomogeneous densities and the fine details along 
the view direction. 

This integration is done by considering piecewise linear segments 
with samples taken on the cloud surface as shown in Figure 10 
(f.e., exponentially spaced samples are read along the view direc¬ 
tion). For one segment z, the single scattering formula gives 

Ti = Mie(0) f 

J Xi 

= ( 1 ) 

where Mie is the Mie phase function and the extinction coeffi¬ 
cient defined in Equation 4, Appendix A. Note that since the Mie 
phase function is wavelength-dependent, we encode and use Mie as 
an RGB function. In multiple scattering, this dependence is blurred 
out, thus T needs only be a scalar. 



Figure 10: Piecewise linear integration of the single scattering paths. 

6.3 Opacity 

The opacity a is computed by integrating the extinction along the 
view direction. This results 'ma — e~'^^dx. Since the extinction 

coefficient k is varying in space, we discretize this integration by 
sampling segments along the view direction as for single scattering. 
For each segment z, we read k in the cloud volume and accumulate 
the opacity 



This is a standard ray marching procedure. 

7 Implementation on the GPU 

At rendering time, the cloud mesh is rendered using a fragment 
shader that computes the light transport from the lit surface as ex¬ 
plained in Section 6. Since this process is fairly involved we rely 
heavily on GPU capabilities to make it efficient, especially for com¬ 
puting distances to the filtered cloud surface, which we explain in 
this section. 

The 22 intermediate functions A through X described in Sec¬ 
tion 5.2 necessary to compute the canonical transport function are 
discretized and stored in textures, as well as the Adze phase func¬ 
tion used for the single scattering. Overall this gives a few textures 
occupying less than 2MB in video memory. 

To manage efficiently the surface processing (e.g., filtering, com¬ 
puting distances) on the GPU we rely on depth maps. We create 
two depth maps ZmiriL, ZmaxL and a normal map Nl from the 
light point of view. This is done at each frame to allow for the 
animation of the cloud shape and of the light. 

A convenient way to approximate the filtering of the surface ac¬ 
cording to a kernel of size a is to rely on the MIP-mapping of the 
Z values in the depth map. To compute the MIP-map correctly 
(i.e., without taking into account pixels where there is no cloud), 
we add an alpha channel A and we follow the same scheme as for 
an alpha-premultiplied texture. In fact the premultiplication is done 
automatically when setting the default Z value to 0. After the MIP- 
map pyramid (Z°,..., Z^) is computed, a filtered Z value is ob¬ 
tained by as in [Dachsbacher and Stamminger 2003]. 

Computing the distances d (Section 6.1) and k (Equation 1) to the 
filtered surface in the light direction is done by simply accessing 
the associated MIP-mapped depth maps. To compute the thickness 
t, we subtract ZmaxL — ZmiriL- The orientation of the filtered 
surface is obtained by reading the MIP-mapped normal map Nl. 
Thanks to this representation, all the parameters (l)L,dA) 

can be computed and the whole algorithm described in Section 6 
can be implemented in a fragment shader. Eor each set of scattering 
orders the shader first finds the collector (c, a) using our algorithm 
described in Section 6.1, then computes the associated intensity T 
using the compressed canonical transport function described in Sec¬ 
tion 5. 
























8 Volumetric enrichment of the cloud model 

Surface-based cloud approaches generally rely on a transparency 
shader to fake the effect of inhomogeneities on the silhouette. 
Volume-based methods are limited in resolution due to the huge 
memory requirement, which usually makes the silhouette lack de¬ 
tails (see Figure 5). Our approach combines the best of both: we 
define a volumetric Hypertexture in a layer below the surface in or¬ 
der to have high-resolution 3D effects on the cloud edges with low 
memory requirements (see Figure 6). Thus, the volume representa¬ 
tion handles all the frequencies of cloud details that are not handled 
by the cloud mesh. 

The Hypertexture layer is defined through a distance function D{p) 
to the surface (which equals 0 on the surface and increases in¬ 
side the cloud). The procedural details are modulating the density 
iVo (Equation 3) and the extinction coefficient k (Equation 4) with 
y9(p) = S{D{p) + noise{p)) where S' is a sigmoid function and 
noise is a scalar 3D Perlin noise [Perlin 1985]. Due to the cost of 
evaluating distance fields, D{p) is computed and stored in a volu¬ 
metric texture. noiseQ is evaluated on the fiy in the shader. We 
rely on [Crassin and Neyret 2007] to compute and voxelize quickly 
the distance field at high resolution using only the minimal neces¬ 
sary memory usage. 

This volumetric enrichment is used in three parts of the shader. 
When computing the cloud opacity (Section 6.3), we perform a ray 
marching on the GPU [Crassin and Neyret 2007] through this layer 
to integrate the cloud density in this part of the cloud. When com¬ 
puting single scattering (Section 6.2), the extinction coefficient n 
(Equation 4) is also modulated for each segment according to p{p). 
Einally, to account for these details in the multiple scattering com¬ 
putations (Section 6.1), the depth map ZmiuL is modulated: we 
store in ZmiuL the depth of the first voxel having p > po. 

Note that for most pixels except on the silhouette and thin parts, 
strong opacity will stop the ray close to the entry point, so the aver¬ 
age marching length is limited. 

8.1 Final rendering 

The precomputation for the current frame is limited to rebuilding 
the depth and normal maps. This allows us to change the viewpoint, 
lighting conditions, and cloud shape. 

Terrain and opaque objects are drawn first. Then the cloud mesh is 
drawn using our pixel shader through deferred shading. This pixel 
shader: 

- iteratively finds the collectors (c, a) for the 8 sets of scattering 
orders corresponding to sun illumination (Section 6.1) ; 

- computes the light transport T associated to each set (as ex¬ 
plained in Section 5), and multiplies each by the lighting condi¬ 
tions (intensity, color and visibility) ; 

- computes the analytical single scattering (Section 6.2); 

- computes the opacity (Section 6.3). 

We handle environment illumination like sun illumination. A blue 
source is added above the cloud and a brown one below the cloud. 
We use our multiple scattering algorithm with both of these addi¬ 
tional sources. Outdoor scenes require accounting for aerial per¬ 
spective. We rely on the model of [Hoffman and Preetham 2003]. 
Rendering very bright objects such as clouds also require some 
tone mapping management. We use a simplified, unblurred version 
of [Goodnight et al. 2003]. 

Note that our representation allows us to easily account for points 
of view inside the clouds, since our canonical transport function 
accounts for this case. 


9 Results 

Our tests were conducted on a Pentium 4 at 1.86 GHz with a nVidia 
8800 GTS graphics board. All benchmarks were done at resolution 
800 X 600. 

We tested our method against different cloud shapes: a cloud slab, 
an animated stratocumulus layer without procedural noise, and a 
cumulus cloud with procedural noise. The slab allows us to test 
our algorithm and validate directly its results against the canon¬ 
ical transport function. It also allows us to see features such as 
anisotropic refiectance, anisotropic transmittance in thin slabs and 
isotropic transmittance in thick slabs. The stratocumulus layer al¬ 
lows us to test our method against a shape close to the canonical one 
and to validate the handling of animations. It is composed of 130K 
triangles, and the framerate is lOfps. We can see on this exam¬ 
ple all the wanted cloud features: anisotropic scattering, diffusion, 
backscattering, glory, fogbow. The cumulus cloud model allows us 
to test our algorithm on an arbitrary shape. It is composed of 5K 
triangles with a 512^ Hypertexture. We obtain a framerate of 2fps 
on this model. Note that most of the time is spent in evaluating the 
noise on the fly in the hypertexture. Without the procedural noise, 
the framerate rises to lOfps. This model also displays all the sought 
cloud features. 10 iterations of our collector-finding algorithm are 
sufficient in all cases to reach convergence. Eigure 12 and the teaser 
show the result of our method on various cases. 

10 Discussion and Future Work 

We address the limitations of other real-time approaches [Harris 
and Lastra 2001; Riley et al. 2004; Hegeman et al. 2005] by treat¬ 
ing accurately light paths of all orders and by using a more detailed 
cloud shape. We use the Mie phase function for cloud droplets, 
which provides us with the right anisotropy and enables visual ef¬ 
fects such as glory and fogbow. Contrary to [Schpok et al. 2003], 
we account for the procedural details in the light transport calcula¬ 
tion. Although similar details could be attained by other real-time 
approaches [Harris and Lastra 2001; Riley et al. 2004; Hegeman 
et al. 2005] by increasing the volumetric resolution of their models 
(at the expense of speed and within memory limits), our method 
brings features that are not handled by these approaches such as 
accurate backscattering and diffusion. We allow for the animation 
of viewpoint, light direction and cloud shape. In case of cloud ani¬ 
mation, we rely on the algorithm described in [Crassin and Neyret 
2007] to recompute in real-time the distance field used by our pro¬ 
cedural noise. 

Like the other methods using a depth-map to solve scattering prob¬ 
lems [Dachsbacher and Stamminger 2003], this representation has 
some issues regarding resolution and non-convex shapes. The pre¬ 
cision of this representation is limited by the resolution of the depth 
map. As a result, full scalability of the method would probably re¬ 
quire alternative depth maps techniques [Stamminger and Drettakis 
2002; Arvo 2007]. In case of non-convex shapes, light transport is 
slightly underestimated in shadowed regions (see Eigure 11(a)): the 
algorithm incorrectly assumes that there is matter between the shad¬ 
owed surfaces and the lit surface. In consequence, light is consid¬ 
ered more attenuated than it should be. In practice, this error hap¬ 
pens on shadowed parts of the cloud, where environment (Le., sky 
and ground) illumination is predominant. As mentioned in [Dachs¬ 
bacher and Stamminger 2003] who suffer from the same issues, this 
issue can be solved by using a depth peeling technique to treat non- 
convex shapes as a collection of convex shapes. Using deep shadow 
maps [Lokovic and Veach 2000] would also help addressing this is¬ 
sue and increase the accuracy of single scattering. 

One limitation of our approach is the assumption that light arriving 


to the viewer goes through only one connected collector {i.e., one 
most probable path) per set of scattering orders. As shown on fig¬ 
ure 11(b) it is not always the case. In these configurations, our 
algorithm accounts for only one collector, thus underestimates the 
amount of light transmitted. Note that since we look for several 
collectors (one for each of the 8 sets of scattering orders) per pixel, 
this error concerns only a fraction of the pixel color. One solution 
might be to look for several collectors per set of scattering orders 
with different initial values. 

Searching for a collector c corresponds to searching for the entry 
point c on the lit surface of a local maximum of the light transport 
T. Our iterative algorithm corresponds to a fixed-point method: we 
look for the fixed point /(c) = c where / represents our canonical 
transport function. As future work, better techniques from the field 
of optimization could be used to ensure faster convergence and to 
handle the case of several maxima of T. 

Our light transport computation on GPU using our collector-finding 
algorithm {i.e., the main contribution of this paper) is fast enough 
for interactive applications (lOfps). The speed of this method can 
be still be improved. Moreover, the computation of the procedu¬ 
ral noise and the GPU ray marching in our implementation is an 
unoptimized version of [Crassin and Neyret 2007] and is computa¬ 
tionally expensive (80% of the rendering cost). 

As future work, we would like to take into account higher-scale 
light effects such as interrefiections between clouds and cloud 
lobes. Also, the method described in [Bouthors et al. 2006] can 
be applied to this work to compute the interrefiections between the 
clouds and the ground, which have been shown to be of visual im¬ 
portance. 


V 

(a) (b) 

Figure 11: (a): Limitation of using depth maps. The radiance of the 
red surface is slightly underestimated. However, this error is of low visual 
importance, (b): Limitation of searching for only one collector per set of 
scattering orders. In this example two collectors of equal importance exist, 
but our method will find only one. 

11 Conclusion 


resentations {e.g., meshes) can be used for the collector-finding al¬ 
gorithm, which would remove the limitations due to depth maps. It 
could also be used to find starting paths in a volumetric Metropolis 
light transport simulation. 
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A Measuring light transport in a cloud slab 

As explained in Section 2, the appropriate phase function for clouds 
is Mie scattering. It is a complex oscillating function depending 
of droplet size and wavelength. Like [Bouthors et al. 2006], we 
preintegrate the effect of the DSD on the Mie phase function and we 
rely on the modified Mie approximation in which the narrow (5°) 
strong forward peak is suppressed out of the phase function and 
converted into 50% lowered density [Lenoble 1985]. This proved 
to yield equivalent results while needing to consider only half of the 
scattering events and gave a better conditioned phase function. 

For the DSD, we draw on the modified Gamma distribution [Levin 
1958] (adapted to cloud droplets) defined by 



This function describes the density N(r) of droplets of radius r, 
with Tn the characteristic radius of the distribution, 7 represent¬ 
ing its broadness, and A'o is the total density of droplets. F is 
the gamma function, p is a modulation factor for the total den¬ 
sity brought by our procedural enrichment described in Section 8 . 
A cloud droplet size distribution is then only described by the three 
parameters Nq, Vn and 7 . In terms of optical properties, the effec¬ 
tive radius corresponding to this distribution is re = (7 + 2 )rn. 
We consider the following typical parameters for cumulus clouds: 
re = 6/im, ^ — 2, Nq — 4.10®m“^. 

The extinction coefficient n is defined by 

n — pNoTtrl. (4) 




We have proposed a study of light transport yielding a new formu¬ 
lation for the macroscopic behavior of light in participating media 
that is suited for computer graphics. We have proposed a new way 
of computing radiative transfer through a volume of homogeneous 
participating media with a dedicated optional method to handle in¬ 
homogeneous boundaries. Contrary to previous approaches, our 
method only needs to walk the cloud boundaries rather than walking 
through the whole volume. We have demonstrated this technique 
on detailed cumulus-type clouds. As shown on our results, we cor¬ 
rectly reproduce visual features such as back-scattering, anisotropic 
multiple scattering, glory, etc. 

We believe this approach can be used in other applications where 
multiple scattering in well-bounded participating media is impor¬ 
tant such as sub-surface scattering, e.g., as an alternative to the 
dipole approximation [Jensen et al. 2001]. As future research, our 
approach could also be used in offline rendering, where other rep¬ 


It yields the extinction function which is the probability to 
traverse the cloud along a path of length x without hitting a droplet. 
It is used in all scattering equations, including opacity (Equation 2), 
single scattering (Equation 1) and multiple scattering (Equation 5). 

The multiple scattering equation has been extensively discussed in 
the literature [Chandrasekhar 1960; Kajiya and von Herzen 1984]. 
It can be written in the form 

=T(p,i!)+/ [ PiCo-^')Tip,a')duj' (5) 

tv dp 47r 

for the intensity T reaching a point p in direction cJ, with P{0) the 
phase function of the media (here. Me). This integro-differential 
equation can be solved computationally through Monte-Carlo inte¬ 
gration. The contribution of each order of scattering in the intensity 
T, as well as the collector information (c, cr) can be easily tracked 
during the integration. For our light transport characterization step 






(Section 5), we computed the solution of Equation 5 in a plane- 
parallel slab for various values of {(^v, o?, t) by ray tracing 

multiple scattering from the eye. We stored the contribution T of 
each scattering order into a database, along with the collector data 
(c, cr). These computations were done with an error margin of 5% 
at the 95% level. We computed them against 10 million different 
sets of values for o?, t), resulting in a raw database 

size of 25GB. These computations took several weeks on a 100- 
nodes, dual-core 900MHz Itanium-2 cluster. The database and fit¬ 
ting results will be made accessible online. 

The analysis of these results {i.e., finding lower-dimensional func¬ 
tions fitting the results) was done empirically by plotting (T, c, a) 
against the input parameters and looking for remarkable behaviors. 
It was inspired by previous approaches such as Chandrasekhar’s X- 
and Y -functions. The fitting itself was done using classical non¬ 
linear least square optimization with MATLAB. 
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Figure 12: Some results of our method. 














